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Abstract 

This paper considers stochastic hybrid stress quadrilateral finite element anal¬ 
ysis of plane elasticity equations with stochastic Young’s modulus and stochastic 
loads. Firstly, we apply Karhunen-Loeve expansion to stochastic Young’s modulus 
and stochastic loads so as to turn the original problem into a system containing a 
finite number of deterministic parameters. Then we deal with the stochastic field and 
the space field by fc—version/p—version finite element methods and a hybrid stress 
quadrilateral finite element method, respectively. We show that the derived a priori 
error estimates are uniform with respect to the Lame constant A 6 (0,+oo). Finally, 
we provide some numerical results. 

Keywords, stochastic plane elasticity Karhunen-Loeve expansion hybrid 
stress finite element k x h —version p x h —version uniform error estimate 


1 Introduction 

Let D C R 2 be a bounded, connected, convex and open set with boundary dD = 
8Dq U 8D\ and meas(3Z?o) > 0, and let {£l,F,V) be a complete probability space, where 
P, J 7 , V denote respectively the set of outcomes, the cr-algebra of subsets of P and the 
probability measure. Consider the following stochastic plane elasticity equations: for 
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almost everywhere (a.e.) 9 £ Q 


—diver(-, 6) = f(-, 9), in D, 

< cr(-,9) = Ce(u(-,0)), in D, (1-1) 

u (' 5 9)\dD Q = 0, <x(-, 9)n\g Dl = g(-, 9), 

where a : D x 0 — > denotes the symmetric stress tensor field, u : D x Q — > R 2 the 

displacement field, e(u) = (yu + v T u)/2 the strain with y = (gfy ~&^) T for x = {xi,x 2 ), 
f : D x —>• i ? 2 the body loading density and g : dD\ x Q —> R 2 the surface traction, n 
the unit outward vector normal to dD, C the elasticity modulus tensor with 


Ce(u) = 2^e(u) + Adivul, 


I the 2x2 identity tensor, and p, A the Lame parameters given by fj, = 2 (u-u) > ^ = 
(l+ufii-^u) f° r P^ ane strain problems and by /./ = A = ^ or pl &ne stress 

problems, with v £ (0, 0.5) the Poisson ratio and E : D x 12 —» R the Young’s modulus 
which is stochastic with 


0 <C G-min — ^) — G max a.e. in D x (1.2) 

for positve constants e m j n and e rnax . Since in the analysis of this paper we need to use an 
explicit form of E, we rewrite the second equation of (1.1) as 

a(-,9) = ECe(u(;9)), (1.3) 

where the tensor C := iC depends only on the Poisson ratio u. 

It is well-known that the standard 4-node displacement quadrilateral element (abbr. 
bilinear element) yields poor results for deterministic plane elasticity equations with bend¬ 
ing and, for deterministic plane strain problems, at the nearly incompressible limit. To 
improve its performance, Wilson et al. [26, 24] developed methods of incompatible modes 
by enriching the standard (compatible) displacement modes with internal incompatible 
displacements. Pian and Sumihara [17] proposed a hybrid stress quadrilateral element 
(PS element) based on Hellinger-Reissner variational principle, where the displacement 
vector is approximated by isoparametric bilinear interpolations, and the stress tensor by 
a piecewise-independent 5-parameter mode. Xie and Zhou [31, 32] derived robust 4-node 
hybrid stress quadrilateral elements by optimizing stress modes with a so-called energy- 
compatibility condition, i.e. the assumed stress terms are orthogonal to the enhanced 
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strains caused by Wilson bubble displacements. In [35] Zhou and Xie gave a unified anal¬ 
ysis for some hybrid stress/strain quadrilateral methods, but the upper bound in the error 
estimate is not uniform with respect to the Lame parameter A. Yu, Xie and Carstensen 
[33] derived uniform convergence results for the hybrid stress methods in [17] and [31], in 
the sense that the error bound is independent of A . 

In the numerical analysis of stochastic partial differential equations, stochastic finite 
element methods, which employ finite elements in the space domain, have gained much 
attention in the past two decades. In the probability domain, the stochastic finite element 
methods use two types of approximation methods, statistical approximation and non- 
statistical approximation. Monte Carlo sampling(MCs) is one of the most commonly used 
statistical approximation methods [22], In MCs, one generates realizations of stochastic 
terms so as to make the problem deterministic, and only needs to compute the determin¬ 
istic problem repeatedly, and collect an ensemble of solutions, through which statistical 
information, such as mean and variance, can be obtained. The disadvantage of MCs lies 
in the need of a large amount of calculations and its low convergence rate. There are 
also some variants of MCs such as quasi Monte Carlo[6] and the stochastic collocation 
method[2, 14, 15, 16]. 

Non-statistical approximation methods mainly contain perturbation methods, Neu¬ 
mann series expansion methods[10] and so on at the beginning. But these methods are 
limited to the magnitude of uncertainties of stochastic terms and the accuracy of calcu¬ 
lation. Later, polynomial approximation is used for the stochastic part. For example, 
Polynomial chaos (PC) expansion is applied in [27, 10] to represent solutions formally 
and obtain solutions by solving the expansion coefficients [9, 13]. Generalized polynomial 
chaos (gPC) is used to express solutions in [12, 28, 29]. According to [30], one can achieve 
exponential convergence when optimum gPC is chosen. Subsequently, it was further gen¬ 
eralized [1, 7] that p version, k version and p-k-version finite element methods could be 
used for the approximation of the stochastic part. 

So far, there are very limited studies on the numerical solution of the stochastic plane 
elasticity equations (1.1). In [11] a generalized nth order stochastic perturbation technique 
is implemented in conjunction with linear finite elements to model a ID linear elastostatic 
problem with a single random variable. In [9] the numerical solution of problem (1.1) 
is considered with stochastic Young’s modulus E, where PC approximation and bilinear 
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finite elements are applied respectively to the stochastic domain and the space domain. 
We refer to [5, 25] for some other related studies. In this contribution, we shall propose and 
analyze stochastic k x h —version and p x h —version finite element methods for the problem 
(1.1), where we use k— version/p—'version finite element methods for the stochastic domain 
and PS hybrid stress quadrilateral finite element for the space domain. 

We arrange the paper as follows. In Section 2 we show stochastic mixed variational 
formulations of (1.1), and give the existence and uniqueness of the weak solution. Section 
3 discusses the approximation of the stochastic coefficient and stochastic loads, as well as 
the truncated stochastic mixed variational formulations. Section 4 analyzes the proposed 
stochastic k x h —version and p x h —version finite element methods and derives uniform a 
priori error estimates. Finally, Section 5 provides some numerical results. 

2 Stochastic mixed variational formulations 

2.1 Notations 

For the probability space (fi, F, V) and an integer m, denote 
Lp(0) := |V| Y is a random variable in (0, F,V) with J \Y(6)\ m dP(9) < +oo 
If Y € Lp(Q), we denote its expected value by 

E[Y} = f Y{6)dP{6) = f ydF(y), (2.1) 

Jn Jr 

where F is the distribution probability measure of Y, given by F(B) = P(Y~ 1 (B)) for any 
borel set B in R. Assume that F(B ) is absolutely continuous with respect to Lebesgue 
measure, then there exists a density function for Y, p : R — > [0, +oo), such that 

E[Y} = [ yp{y)dy. (2.2) 

Jr 

We denote by H m (D) the usual Sobolev space consisting of functions defined on the 
domain D, with all derivatives of order up to m square-integrable. Let (■,-)i? m (n)be the 
usual inner product on H m (D). The norm || • || m on H m (D ) deduced by (•,-)i^ m (D) is 
given by 

I Him := ( ^ Mj) 1 / 2 with the semi-norm \v\j := ( ^ ||H a u||Q) 1//2 . 
o <j<rn |a[ =j 
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In particular, L 2 (D) := H°(D). Denote 

L°°(D ) := {re : HreHoo := esssup xe D\w(x)\ < oo}. 

We define the following stochastic Sobolev spaces: 

L 2 p(Q-,H m (D)) := {re : re is strongly measurable with w(-,0) € H m {D) for 8 € D and ||re||, 7 i < + 00 } 

Lp(D;L 00 (D)) := {re : re is strongly measurable with w(-,9) € L°°(D ) for 8 € fi and ||re ||55 < + 00 } 
where the norms || • ||m> || • ||oo are respectively defined as 

IMIm : = IMloo : = esssu P6en\\ w ('i 0)||oo. (2.3) 

On the other hand, since stochastic functions intrinsically have different structures with 
respect to 9 € D and x G D, we follow [1] to introduce tensor spaces for the analysis of 
numerical approximation. Let -Xi(fi), W)(-D) be Hilbert spaces. The tensor spaces X\ (D)O 
X 2 (D) is the completion of formal sums < f>(d,x) = £fc=i t ... tn Ui(8)vi(x),Ui € Xi(Q),Vi € 

X 2 (D), with respect to the inner product^,</>)xi®x 2 := '^i,j(ui,Uj)x 1 (vi,Vj)x 2 - Then, 
for the tensor space L 2 p ( D) tg> H m (D), we have the following isomorphism: 

L 2 P (Sl- H m {D )) ~ L 2 p (Q) 0 iL m (D). 

For convenience, we use the notation a < b to represent that there exists a generic 
positive constant C such that a < Cb, where C is independent of the Lame constant A 
and the mesh parameters h, k, the polynomial degree p in the stochastic k x h— version 
and p x h— version finite element methods. 


2.2 Weak formulations 


Introduce the spaces 


: = 


L 2 (D; R 2 S y%) := {t : D ->■ R 2x2 \ T zj € L 2 (D), nj = Tp, i,j = 1, 2}, if meas(<TDi) > 0, 


V D :={veH\D) 2 :v\ dDo = 0} 

4™) I? 2x2 

{t € L 2 (D; R 2 ym) : / D trrdx = 0 with trace trr := rn + T 22 }, if <9Z?i = 0. 
Then the weak problem for the model (1.1) reads as: Find (<x, u) €E L 2 P (V,] Ep) x 
L 2 p (p,\ Vo) such that 

a(cr, r) — 6(r, u) = 0, Vr € L 2 P (Q] Ep>), 

6(tr, v) = £(v), Vv <E L 2 p (n;V D ), 


(2.4) 
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where the bilinear forms a(-,-) : L P (Q; Ed) x Lp(H; Ed) —> R, &(■,•) : Pp(I2; Ed) x 
L p(fl; Vd) —> R and the linear form l : Lp(H; Vd) —> P are defined respectively by 

a(er,r) := E[ f — a : C _1 rdx] = / f — cr : C _1 rdxdP(0), (2-5) 

Jd P Jd Jd P 

b(r,u) := E[ f r : e(u)dx] = f f r : e(u)dxdP(0), (2-6) 

Jd Jo Jd 

£(v) := E[ f fvdx + f g • vds] = f f fvdxdP(0) + f f g-vdsdP(0). (2.7) 

ifl JoD x JO JD JO JdD\ 

Here a :t = £- j=1 cr^T^.. 

It is easy to see that the following continuity conditions hold: for er,r G Tp(fl; Ed), 
vGL^H; Vd), 


a(o-,T) < ||<t|| 5 ||t||q, 6(t,v) < ||r|| 5 |v| T , £(v) < (||% + llgllo,^) Ml- ( 2 -8) 

According to the theory of mixed finite element methods [3] [4], we need the following two 
stability conditions for the well-posedness of the weak problem (2.4): 

(A) Kernel-coercivity: for any r € Z° := {r G Lp(S2; Ed) : b(r, v) = 0, V v G 
Lp(H; Vd)} it holds 

IMI§ ^ (2.9) 


(B) Inf-sup condition: for any v G L p (Q; Vd) it holds 


|v|i < 


sup 

0^TGLp(£2; Sc) 


b(r , v) 


Theorem 2.1. The uniform stability conditions (A) and ( B) hold. 


( 2 . 10 ) 


Proof. For any r G Z°, we have, a.e. 6 G Q, r(-, 0) G (r G Ed : /d T ' e ( v )dx = 0 V v G 
Vd}. According to Theorem 2.1 in [33] and the assumption (1.2), it holds 

/ t(-, 0) : r(-,0)dx < f -Lr(-, 0) : C _1 (-, 0)r(-, 0)dx, 

Jd Jd P 

which leads to 


O JD 


r : rdxdP(0) < [ [ i • r : C " 1 rdxdP(0), 
Jo Jd P 


i.e. (A) holds. 

Let v G Lp(H; Vd) and notice e(v) G Tp(fl; Ed)- Then 

, , m ^ Jo Id t : e(v)dxdP(0) 

k(v)lo< sup -n—-• 

tgl 2 p (q-, s d )\{o} ll T llo 

Hence (B) follows from the equivalence between the two norms |e(v )|q and |v|y on Lp(fi; Vd). 

□ 
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In view of the above conditions, we immediately obtain the following well-posedness 
result: 

Theorem 2.2. Assume that f G L 2 p (Q, L 2 (D) 2 ), g G L 2 p (TL,L 2 (dD\) 2 ). Then the weak 
problem (2-4) admits a unique solution ( cr,u ) G Lp(f2; £d) x L 2 p (Q\ Vo) such that 

IMIo + Mi ~ ll/llo + Iloilo, dD 1 - (2- 11 ) 


3 Truncated stochastic mixed variational formulations 

In order to solve the weak problem (2.4) by deterministic numerical methods, we firstly 
approximate the stochastic coefficient E and the loads f, g by using a finite number of 
random variables; we refer to [21] for several approximation approaches. Here, we only 
consider the Karhunen-Loeve(K-L) expansion. 


3.1 Karhunen-Loeve(K-L) expansion 


For any stochastic process 0(x, 6) G L 2 p {Q ,; L 2 (D)) with covariance function cou[0](xi, X 2 ) : 
D x D —^ R , which is bounded, symmetric and positive definitely. Let {(A n , 6 n )}^ =1 be 
the sequence of eigenpairs satisfying 


/ cov [<j)\ (xi,x 2 ) 6„(x 2 ) dx 2 = A n 6 n (xi), 
J D 


(3.1) 


+00 „ 

VA n = / cou[0](x,x)dx, 
n=1 ■ /D 


/D 


6 i(x) 6 j(x) dx = i,j = 1,2,--- , (3.2) 

and Ai > A 2 > • • • > 0. Then the Karhunen-Loeve(K-L) expansion of </>(x, 0) is given by 

OO 

</>(x, 0) = -E[<£](x) + ^ \/A nb n (x.)Y n (6), (3.3) 

n=l 

and the truncated K-L expansion of </>(x, 6) is 

N 

</>Ar(x,0) = £[<?!>](x) + ^ V / A nb n (-x)Y n (0). (3.4) 


n= 1 


Here {y r) ,}^ =1 are mutually uncorrelated with mean zeros and unit variance with Y n (6) = 
7 ^ / D (^)(x, 0) - E[(/)] (x)) 6 n (x)dx. 

By Mercer’s theorem [20], it holds 

+OO 

sup E[(cp — 0tv) 2 ](x) = sup A n fe 2 (x) —> 0. as IV —> oo. (3.5) 


xg.d 


xG-D 


n=7V+l 
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In what follows we show the estimation of the truncated error cj) — in norms 11 ■ | |g 
and || • |jog, respectively. 

From (3.2) it follows 

+OO 

||0-<M|~= ^ n and 11^ - 0iv|lo -> 0 as iV->+ OO. (3.6) 

n=N+l 

Obviously the convergence rate of \\<f> — 4>n\\q is strongly depending on the decay rate of 
the eigenvalues X n , which ultimately depends on the regularity of the covariance function 
cov[(/)\. Generally, the smoother the covariance is, the faster the eigenvalues decay, which 
implies the faster \\<f) — </>jv||q converges to zero. Now we quote from [23] the following 
definition (Definition 3.1, which are related to the regularity of cov[(j)]) and lemma (Lemma 
3.1, which gives the decay rate of the eigenvalues A n ). 

Definition 3.1. [23] The covariance function cov[4>\ : D x D —>■ R is said to be piecewise 
analytic/smooth on DxD if there exists a finite family (Dj)i<j<j C R 2 of open hypercubes 
such that D C U j = 1 Dj, Dj n Dy = 0, Vj 7^ j' and cov[c[>] \ o. t x d : i , has an analytic/smooth 
continuation in a neighbourhood of Dj x Dji for any pair 

Lemma 3.1. [23] If cov[(j)\ is piecewise analytic on DxD, then for the eigenvalue sequence 
{A n }n>i, there exist constants c,\ , depending only on cov[4>] such that 

0 < A n < cie _C2nl/2 , Vn > 1. (3.7) 

If cov[(j)\ is piecewise smooth on DxD, then for any constant s > 0 there exists a constant 
c s depending only on cov[4>] and s, such that 

0 < \ n < c s n~ s , Vn > 1. (3.8) 

By Lemma 3.1, we immediately have the following convergence results. 

Lemma 3.2. If cov[4>\ is piecewise analytic on D x D, then there exists constants C\,C 2 
depending only on cov[4>\ such that 

\\<[-cf N \\ d <^(l + c 2 N 1 / 2 )e-^ Nl/2 , \/N> 1. (3.9) 

C 2 

If cov[(j)] is piecewise smooth on D x D, then for any s > 0 there exists C s depending only 
on cov[4>\ and s, such that 


H0-<Mlo < C S N~ S 


VN > 1. 


(3.10) 



To estimate || <j) — 4>n\\^5, we make the following assumption: 

Assumption 3.1. The random variables {Y n (0)y^ =1 in the K-L expansion are indepen¬ 
dent and uniformly bounded with 

\\Yn(0)\\L°°(n)<C Y , Vn > 1, 

where Cy is a positive constant. 

Lemma 3.3. [8, 23] Suppose Assumption 3.1 holds. If cov[(j)] is piecewise analytic on 
D x D, then there exist a constant c > 0 such that, for any s > 0, it holds 

\\<f> - 4>n\ loo < Ce- c ^ 2 ~ s ^ Nl/2 yN > 1, (3.11) 

where C is a positive constant depending on s,c , cov[(j)\ and J given in Definition 3.1. If 
cov[f)\ is piecewise smooth on D x D, then for any t > 0, r > 0, it holds 

I \<t> - </W||oo < C'N^- T V\\/N > 1, (3.12) 

where C’ is a positive constant depending on t,r, cov[(j>\ and J. 

Remark 3.1. We note that we need to solve the integral equation (3.1) to obtain the 
K-L expansion (3.3). For some special covariance functions, the equation can be solved 
analytically [10], but for more general cases numerical methods are required [8, 18, 23]. 


3.2 Finite dimensional approximations of E, f, g 


In this section, we use the K-L expansion to approximate E, f and g. 
For E, assume its truncated K-L expansion is of the form 


N 


E N (x,0) = E N (x,Y 1 {6),...,Y N (d)) = E[E\{x) + 


A n b n (x)Y n (6), 


(3.13) 


n= 1 


where {(A,,., b n (x))}]] =1 and {Y n (6)}]] =1 are the corresponding eigenpairs and random vari¬ 
ables, respectively. 

As for f = (/i, f 2 ) T and g = (gi,g 2 ) T , we need to apply the K-L expansion to each of 
their components. In this paper, following similar ways as in [1, 2] to avoid use of more 
notations, we assume the truncated K-L expansions of f and g take the following forms: 


fv(x, 0) = f]v(x, Y\[Of ..., Kv(0)) = 


flN 

flN 


ElfiK*-) 

£[/ 2 ](x) 


JV 

E 

n =1 


Aln^ln(x) 

A 2 n6 2 n(x) 


Y n (0), 

(3.14) 
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g w (x, 9) = g Af (x,Fi(0), Y n {6)) 


9\n 

g2N 


E [gi] ( x ) 

-% 2 ](x) 


N 


+ £ 


) Y„(0), 
vA 2 n 62 „(x) y 

(3.15) 


where {(A in ,6j n (x))}^ =1 , {(A in ,M x))}^ =1 ,* = 1,2 are the corresponding eigenpairs. 


Remark 3.2. /n practice, the Young's modulus E, the body force f and the surface load g 
may be independent. In such cases, the random variables {Y n (6)}^ =1 in the truncated K-L 
expansions (3.13)-(3.15) for E, f\, f 2 , g\,g -2 may be different from each other. However, 
the analysis of this paper still applies to these cases. 


3.3 Truncated mixed formulations 


By replacing E, f, g with their truncated forms Ejy, f/v, gj\r hi the bilinear form af, •), 
given in (2.5), and the linear form £(■), given in (2.7), we can obtain the following modified 
mixed variational formulations for the weak problem (2.4): find (<xat,Utv) € Lp(12; Ep) x 
Lp(fi; Vo) such that 


a N (a N ,r) - b(r,u N ) = 0, Vr G L 2 p (tl; Sp), 
b(a N , v) = ^jv(v), Vv G Lp(12; Vp). 


(3.16) 


We recall that {Y) l (0)}^L 1 are the random variables used in the K-L expansions of E, 
f and g, which are assumed to satisfy Assumption 3.1. In what follows we denote 

N 

Y := (n,^ s ,..,^v), T n := K„(f2) C R, r := r„, (3.17) 

n=l 


and let p : T —> R be the joint probability density function of random vector Y with 
p G L°°(r). According to Doob-Dynkin lemma [19], the weak solution of the modified 
problem (3.16) can be described by the random vector Y as 


Uat(x, 6) = Uat(x, Y), crjyfx., 9) = ct n (x, Y), 


and, by denoting y := (yi,y 2 , • • • , j/jv), the corresponding strong formulation for (3.16) is 
of the form 


—diver tv(x, y) = fjv(x, y), V(x, y) G D x T, 

er A r(x, y) = E N C e(ujv(x, y)), V(x, y) G D x T, 

< 

ujv(x, y) = 0, V(x, y) G 8D 0 x F, 

er A r(x, y)n = g JV (x, y), V(x, y) G ODi x T. 
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Recall that p : T — y R is the joint probability density function of random vector Y . 
We introduce the weighted L 2 -space 

L 2 p {T) := {v : T —> R j J pv 2 dy < +oo}. (3.19) 

We note that from the norm definition (2.3) it follows 


IMlI = J^p(y)\\ w (-,y)\\m d y = Vu; € L 2 p (T)®H m (D). (3.20) 

It is easy to see that the modified problem (3.16) is equivalent to the following deterministic 
variational problem: find (ctn, un) £ (L 2 (T) <g> Yd) x (L 2 (T) <g> Vo) such that 


a N {a N , r) - b N (r , u p ) = 0, Vr € L 2 (V) ® S D , 
b N (a N , v) = £ n (v), Vv G L 2 (T) <g> V D , 


(3.21) 


where 

a N {a N ,r) := [ p( y) [ J- • cr N : C _1 rdxdy, (3.22) 

Jr Jd E n 

b N {r , utv) := f p{ y) f r : e(uAr)dxdy, (3.23) 

ir Id 

^iv(v) := f p( y) I f/vvdxdy + f p{ y) f g N ■ vdsdy. (3.24) 

Jr Jd Jr JdD 1 

The significance of the form (3.21) lies in that it turns the original formulation (2.4) 
into a deterministic one with perturbations of the Young’s modulus E. the body force f 
and the surface load g. Lemma 3.4 shows, if the perturbations or the truncated errors are 
small enough, we can numerically solve the deterministic problem (3.21) so as to obtain 
an approximate solution of the original problem (2.4). 


Remark 3.3. In some applications it may be more efficient to numerically solve the 
problem (3.21) just in a subdomain TcT, as, of course, will cause that the corresponding 
approximation solution has no value in T\T. 

Lemma 3.4. Suppose that Assumption 3.1 holds and the covariance function, cov[E], of 
E is piecewise smooth (cf. Definition 3.1). Then, for sufficiently large N, the modified 
weak problem (3.16), or its equivalent problem (3.21), admits a unique solution ( ajy , un) £ 
(L 2 (r) <g> Ejj) x (L 2 (r) <g > Vd) such that 

Ik - °w|lo +1«- m\i < We-EnW^ ■ Ikllo + II/-/jvIIo + Ik - ftvllo.aiv ( 3 - 25 ) 

where ( cr,u ) £ L 2 P {D\ S/j) x L 2 p {fl\ Vd) is the solution of the weak problem (2-4). 
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Moreover, (i) if the covariance functions cov[E\, cov[f\ and cov[g\ are piecewise ana¬ 
lytic, then there exists a constant r > 0, and a constant C r > 0 depending only on cov[E], 
cov[f\, cov[g\ and r, such that 

||cr - crjvlI q + \ u ~ u n\j < C r N l / 2 e~ rNl/2 . (3.26) 

(ii) If cov[f\ and cov[g\ are piecewise smooth, then for any s > 0, there exists C s > 0 
depending only on cov[E], cov[f\, cov[g\ and s, such that 


||cr - crjvl |q + \u-u N \j < C S N s . (3.27) 

Proof. We first show the modified problem (3.16) is well-posed. Since the uniform stability 
conditions for the bilinear form &(-,-) and the linear form In{') hold, it suffices to show 
that E]\f is, for sufficiently large N, uniformly bounded with lower bound away from zero 
a.e. in D x ft. In view of Lemma 3.3 and the assumption (1.2), there exists a positive 
integer Nq such that, for any N > Nq, it holds 

e' min <E N < e' max a.e. in D x to, (3.28) 

where e' min and e' max are two positive constants depending only on the bounds of E, i.e. 
e mm and e max in (1.2). Thus, the corresponding uniform stability conditions of the bilinear 
form ajv( - , •) follow from those of a(-,-). As a result, the weak problem (3.16) admits a 
unique solution (<x/v,u/v) € Lp(fl; Eo) x Lp(fl; Vo) with the stability result 


°W|lo + I u ^v|t ^ I|fjv|In + llgjvll 


OydDl 


(3.29) 


for N > Mo- 

Next we turn to derive the estimate (3.25). Subtracting the corresponding equations 
in (2.4) and (3.16), we have 


a N (a - <tjv, r) - b(r, u - uat) = a N (a, r) - a(a, r), Vr € L 2 p (Q; S D ), 
b(a - cr N , v) = £(v) - i N (v), Vv € V D ). 


(3.30) 


Then the desired estimate (3.25) follows from the corresponding stability conditions. 

By Lemmas 3.2-3.3 and the estimate (3.25), we immediately obtain the estimates 
(3.26)-(3.27). □ 
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4 Stochastic hybrid stress finite element methods 


In this section, we shall consider two types of stochastic finite element methods for the 
truncated deterministic variational problem (3.21): k x h version and p x h version. We 
use the PS hybrid stress quadrilateral finite element [17] to discretize the space field and 
A;—version/p—version finite elements to discretize the stochastic held. 

For convenience we assume that the spacial held D is a convex polygon and the stochas¬ 
tic hied T = Y\n=i r n is bounded (cf. Assumption 3.1). 

4.1 Hybrid stress finite element spaces on the spatial field 

Let Th be a partition of D by conventional quadrilaterals with the mesh size h := 
maxTeTh^T, where Iit is the diameter of quadrilateral T £ Th- Let Aj( x±\x^), 1 < i < 4, 
be the four vertices of T, and T* the sub-triangle of T with vertices Aj_i, Aj, Aj + i (the 
index of A , is modulo 4). We assume that the partition Th satisfies the following ’’shape- 
regularity” hypothesis : there exist a constant £ > 2 independent of h such that, for all 
T £ Th, it holds 

h'T ^ Cpt, (4-1) 

where px '■= mmi<j <4 {diameter of circle inscribed in T ,}. 


^3 



Ai | A 2 I 

Figure 1: The mapping Ft 


Let T = [—1,1] x [—1,1] be the reference square with vertices A t , 1 < i < 4(Fig.l). Then 
exists a unique invertible mapping Ft that maps T onto T with Ft ( A, ) = A^, 1 < i < 4. 
The isoparametric bilinear mapping ( x±,X 2 ) = Fr(x 1 , 2 : 2 ) is given by 

xi = a 0 + a\xi + a 2 xix 2 + 03 ^ 2 , ^2 = ^0 + ^ 1^1 + b- 2 xiX2 + 63 ^ 2 , (4.2) 
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where xi,x 2 € [—1,1] are the local isoparametric coordinates, and 


/ 


a o bo 


\ 


«i h 
a 2 b- 2 
\ as &3 / 


/ i 


1 1 1 




- 111-1 
1 - 11-1 
- 1-11 1 


/ \ 


X 


(2) J2) 


1 


Xn 


X 


(3) „(3) 


1 


Xn 




X 


(4) 14) 


2 / 


In Pian-Sumiharas hybrid stress hnite element (abbr. PS element) method for de¬ 
terministic plane elasticity problems, the piecewise isoparametric bilinear interpolation is 
used for the displacement approximation , namely the displacement approximation space 
Von C Vo is chosen as 

Vdh '■= { v e Vd ■ v = v\ t oF t € span{l,xi,x 2 ,xix 2 } 2 , V T € %}■ (4.3) 

In other words ,for v = ( v,uj) t G Vy h with nodal values v(Tj) = ( Vi,Ui) T on T, v is of the 


form 


where 


v = 


VO + VyX! + V 2 X\X 2 + V 3 X 2 

W 0 + W\xi + W 2 X\x 2 + W 3 X 2 


1 Vo Wo ^ 


( 1 1 1 1 ) 


^ V\ ^ 

Pi Wi 

1 

-1 1 1-1 


V 2 UJ 2 

v 2 w 2 

- 4 

1 -11-1 


v 3 u 3 

\ P 3 W 3 J 


l - 1 - 1 1 1 J 


\ V 4 w 4 y 


To describe the stress approximation of PS element, we abbreviate the symmetric 
Til T12 


tensor r = \ ~ ~ J to r = (rn, t 2 2, ti 2 ) t 

ti 2 t 22 

element takes the following form on T: 


The 5-parameter stress mode of PS 





/ 1 0 0 

x 2 

“39, \ 



T = 

T 22 

= 

0 1 0 


Xl 

(3 T for P T = {Pl,...,PZ) T (E R 5 . 

(4.4) 


v ri2 y 


\ 0 0 1 


°'3 

ai X2 b 3 


xi y 


Then the corresponding stress approximation space for the PS hnite element is 

So/j := {r € T10 : r = t|toTt is of form (4.4), VT € Th}- (4-5) 

4.2 Stochastic hybrid stress finite element method: k x /^-version 

This subsection is devoted to the stability and a priori error analysis for the k x h- 
version stochastic hybrid stress hnite element method ( k x h-SHSFEM). 
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4.2.1 k x /r-SHSFEM scheme 


We first use the same notations as in [1] to introduce a ^-version tensor product finite 
element space on the stochastic field r = n;T=i r n c r n . 

Consider a partition of T consisting of a finite number of disjoint R w -boxes, 7 = 
n^Li(an, tt) with (aZ,bZ) C T n and the mesh parameter k n := max 7 \bZ — aZ | for n = 
1,2, - - - ,N. 

Let q = (qi,q- 2 , ■■■,qN ) be a nonnegative integer muti-index. We define the A:—version 
tensor product finite element space Ij ^ 1 as 

y k : = ®n=i Y k^ Y kn : = {^ : r n ^ R : v\(al,bl) € span{yZ : a = 0, 1, ...,q n ], V 7 } . 

(4.6) 

The k x h-SHSFEM scheme for the original weak problem (2.4), or the modified weak 
problem (3.21), reads as: find ( 07 ^, 11 ^) G (l ^ 1 <g> E r> h ) x (Ik* ® Y Dh) such that 

{ Q'N^G'kh, 1 T~kh) bjy^Tkh, u kh) — 0, VTfc/j G ® ^ 

bN(<Tkh, Vfcft) = ^Af(vfcft), VVfcfc € F^ 1 <g> Vdh- 

Here we recall that 

y k ® E Dfc = spcm{^(y)r(x) : ip G Y^ 1 , r G S Dfc }, 
y k q ® Vbft = span{ip( y)v(x) : G Y k q , v G Vb/J, 
and V Dh , E D h are defined in (4.3), (4.5), respectively. 

4.2.2 Stability 

To show the k x h-SHSFEM scheme (4.7) admits a unique solution, we need some 
stability conditions. We note that the continuity of ajv(-, •), &jv(‘, •) and Cv(') follows from 
their definitions. Then, according to the theory of mixed methods [3], it suffices to prove 
the following two discrete versions of the stability conditions. 

(A h ) Discrete Kernel-coercivity : for any r kh G Z% h := {r k h € Y£®Y, Dh : b N (T k h, v k h ) = 
0, Vv fch € T k q <g> Voh} , it holds: 

lkfc/i||§ ^ a N (T k h,Tkh)- (4.8) 

(B/j) Discrete inf-sup condition : for any G Y ^ 1 tg> Vph , it holds 

I | ^ bjy^TkhiVkh) 

sup —77 - rr —• (4.9) 

o ¥=T kh eY«®X Dh \\ T kh\\ 0 
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To prove the stability condition (A^), we need the following lemma [33]: 


Lemma 4.1. Assume that for any piecewise constant function w, i.e. w £ L 2 (D) with 
w\t = const, VT € Th, there exists v £ Voh with 

2 

0 - 

Then, for any r h £ {r h £ T, Dh : f D r h : e(v h )dx = 0, Mv h £ V Dh }, it holds 

\\Th\\l< [ -i-r h : C-'rhdx. 

Jd E n 

We note that the assumption of this lemma, which was first used in [34] in the analysis 
of several quadrilateral nonconforming elements for incompressible elasticity, requires that 
the quadrilateral mesh is stable for the Stokes element Q1-P0. As we know, the only 
unstable case for Q1-P0 is the checkerboard mode. Thereupon, any quadrilateral mesh 
subdivision of D which breaks the checkerboard mode is sufficient for the uniform stability 

(A,). 



Lemma 4.2. Under the same condition as in Lemma f.l, the uniform discrete kernel- 
coercivity condition ( Aif) holds. 

Proof. For any r k h G Z® h , due to the definitions of spaces Y^ 1 <8> Yi£> h and (g> Vdh we 
easily have T kh (-, y') G \r h £ Y Dh : j D r h : e(v /) )dx = 0, Vv h € V Dh } for any y' £ T. 
From Lemma 4.1 it follows 



T~kh{'-, yOdx < [ — -T kh (-,y'):C 1 T kh (-, y')dx, Vy' £ F, (4.10) 

Jd E n (■ ,y ) 


which immediately implies (A h). 


□ 


To prove the discrete inf-sup condition B/, we need the following lemma: 

Lemma 4.3. For any v kk £ YjJ <g> VDh, there exists r k h G YjJ (8 > such that, for any 

T£T h , 

P(y) f T T kh ■ e{v kh )dxdy= ||r fefe ||| T > ||e(t> fch )||| T . (4.11) 

Proof. The desired result is immediate from Lemma 4.4 in [33]. □ 

Lemma 4.4. The uniform discrete inf-sup condition (B k ) holds. 
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Proof. From Lemma 4.3, for any Vkh £ <g> Vo hi there exists T k h £ l )) 1 < 2 )Y,oh such that 


l|r fc ,||o|v fc ,| T <(E / / T kh ■ T kh dyidy)2 (Y] / p( y) / e(v fc/l ) : e(v fc/l )dxdy) 2 

T Jr Jt t Jr Jt 

E / p(y) / T kh ■ Tfchdxdy < / p(y) / Tkh : e(v fc/l )dxdy, 
p Jr Jt Jr Jr> 

where in the first inequality the equivalence of the seminorm |e(-)|g and the norm || • ||j on 
the space Vd) is used. Then the uniform discrete inf-sup condition (B h) follows 

from 


1 1 ^ ir P(y) It T kh ■ e ( v fc/i)dxdy 

;$-n-iu-< sup 

T' kh &Y*®Y, Dh 


Jr p(y) It T 'kh : e ( v fc/i)dxdy 


'fc/illO 


□ 


In light of Lemma 4.2 and Lemma 4.4, we immediately obtain the following existence 
and uniqueness of the k x h-SHSFEM approximation (a k hi u k h)' 

Theorem 4.1. Under the same condition as in Lemma 4-1 , the discretization problem 
(4-7) admits a unique solution (crkh, Ukh) £ (Y k <8> S^/,) x (Y£ <g) Voh)- 


4.2.3 Uniform error estimation 

In what follows we shall derive a priori estimates of the errors ||er — cr/ c /, | |g and |u — 
u A:/ili which are uniform with respect to the Lame constant A € (0,+ 00 ), where (<r,u) € 
S/j)) x (Lp(f2; Vo)) is the solution of the weak problem (2.4). 

Let (<tat, u/v) £ (Lp(r)(g)X£)) x (Lp(r)<8> Vo) be the solution of truncated weak problem 
(3.21). By triangle inequality it holds 

II* 7 " cr/chllo — II* 7 " <Uv|lo 4” Ht/v <7kh\\o, (4.12) 

|u - U kh \j < |u - u N \j + |uat - Ufcftli, (4.13) 

where the perturbation errors, ||<r — cr^v ||q an d l u — u aHi> are estimated by Lemma 3.4. 
For the finite element approximation error terms ||<r/v — &kh\\Q and |utv — u^p, from the 
stability (A/,), (B/J and the standard theory of mixed finite element methods [3] it follows 

||o , JV-o , fc/i||o + l u JV-u fe / l |r < inf \\<r N -T k h\\o+ inf |ujv-v fc/l L. (4.14) 

r kh £Y*®Y Dh u v kh eY«®V Dh 
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To further estimate the righthand-side terms of the above inequality, we need some 
regularity of the solution (<t/v,ujv)- In fact, it is well-known that the following regularity 
holds: 


kiv(-,y)||i + ||ujv(-,y)|| 2 < ||fjv(-,y)||o + l|gAr(- > y)llo,au 1 , Vy G T. 


(4.15) 


On the other hand, in view of (3.28) and the truncated K-L expansions (3.13)-(3.15), and 
by taking derivatives with respect to y n in (3.18), standard inductive arguments yield 

lffi + 1 <M-,y)llo |^ + 1 UA r(-,y)|i ^ ro _ x 0 „+ ln|f ( \.. r x,| m w_. rr 

- (q + 1)! - 1 - (q + 1)! -~ ^ n > (l|iA r (- 1 y)llo+||giv(',y)llo,0D 1 +l), Vy G I, 

(4.16) 


where 


1 


Hn '■ —max{ / \/A n ||6 n || L cx , ( D ), y A in || 6 in ||o(i — 1) 2), y Am||^m||o, 9 Di(^ — 1)2)}. (4.17) 


Then, thanks to = ®^ = 1 Y)? n and the regularity (4.15)-(4.16), standard interpolation 
estimation yields 

N 

\<rN-Tkh\\o ^ H\\(T N \' 


inf 

Tkh& r ^®^Dh 


< MI„„I| T + l |8 £ + ‘^ll^(r)«^ 


n =1 


N 


< h + ^2(k n -f n ) qn+1 , 


(4.18) 


n =1 


N 


■ t l I ^ . .^ + 1 \\dy n + u A r llL 2 (r)(g)V D 

lujv-VfchlT h\\u N \\^ +2_,( — ) qn+ - 


v^eyc 1 ®v Dh 


n= 1 


(Qn + 1 )! 


N 


< h + ^2(k n -f n ) qn+1 . 


(4.19) 


n =1 


In light of the estimates (4.14) and (4.18)-(4.19), we immediately obtain the following 
conclusion. 


Theorem 4.2. Let ( cr N , u N ) G (L 2 (T) <g> S D ) x (. L 2 (T ) ® V D ) and (a kh , u kh ) G (Yjf ® 
S Dh) x O^k ® Vbh) the solutions of (3.21) and (4-7), respectively. Then, under the 
same condition as in Lemma 4-1 and for sufficiently large N, it holds 

N 

Ikiv - o-kh \\o + | U N - u kh It <h + ^2(k n 7 n) 9n+1 . (4.20) 

n=l 

Remark 4.1. We notice that the estimate (4.34) is optimal with respect to the mesh 
parameters h and k = [k \, k- 2 r - ■ ■ ,kj\f), but not optimal with respect to the polynomial 
degree q = (<?i, • ■ ■ , qN) since it requires k n -/ n < 1 . 
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The above theorem, together with Lemma 3.4, implies the following a priori error 
estimates for the k x /i-SHSFEM approximation {(Tkh^kh)- 

Theorem 4.3. Let (a, u ) € (Lp(0; Ed)) x (L 2 p (LI; Vd)) and ( a kh , «/,-/,.) G (Y^<S>E Dh ) x 
(Yg® v D h ) be the solutions of (2.f) and (f.7), respectively. Then, under the same condi¬ 
tions as in Theorem f.2, it holds 

N 

Ik - <Tkh\\o + l w ~ u kh\i < iV 1 / 2 e _rAri/2 + h + ^2{k n 7„) 9n+1 (4.21) 

71=1 

/or any r > 0 if the covariance functions of E, f and g are piecewise analytic, and holds 

N 

Ik - cr kh Ho + \u- u kh | T < N~ s + h + ^2(k n 7n)' ? " +1 (4.22) 

71=1 

for any s > 0 if the covariance functions of E, f and g are piecewise smooth. 

Remark 4.2. Here we recall that ” < ” denotes ” < C" with C a positive constant 
independent of X , h , N , k. 

4.3 Stochastic hybrid stress finite element approximation: p x h version 

As shown in Section 4.2 and Remark 4.1, the k x /i-SHSFEM is based on the k partition 
of the stochastic field T and requires the mesh parameter k n (n = 1,2, ••• , N) to be 
sufficiently small so as to acquire optimal error estimates. 

In this subsection, we shall introduce a p x h-ve rsion stochastic hybrid stress finite 
element method (p x /r-SHSFEM), which does not require to refine T. We will show this 
method is of exponential rates of convergence with respect to the degrees of the polynomials 
used for approximation. To this end, we first assume 

E n G C°(r,L°°(D))> fjv G C°(T,L 2 (D)), SN €C°(T,L 2 (dD 1 )). (4.23) 

Here 

C°(r,R) := {u : T — > B, v is continuous in y and max||u(y)||s < +oo} (4.24) 

yer 

for any Banach space, B, of functions defined in D. The above assumptions indicate that 
the solution, (tT^v,ujv), of the problem (3.21), satisfies 

a N G C°(F,E d ), u N eC°(r,v D ). 
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Let p := (pi,P 2 , ■■■,Pn ) be a nonnegative integer muti-index. We define the p— version 
tensor product finite element space Z p as 


Z p := < 8 'n=i z n n , := {ip : T n -» R : cp G span{y% : a = 0,1, ...,p n }} 


(4.25) 


(4.26) 


Then the p x /i-SHSFEM scheme reads as: find (cr pft ,u p /j) € (Z p < 8 > S on) x (-^ p ® Vd/i) 
such that 

(.Cphi T~ph) ^N^phi Up/j) — 0, V'T'p/i G Y p 8) E j)j r . 

bN(cr ph , v p/l ) = lN(v p h), Vv,ph G Y p < 8 ) Vb h . 

We note that Z p is a special case of the k —version tensor product finite element 

space Y k q , then, in this sense, the p x /i-SHSFEM can be viewed as a special case of the 

A: x /i-SHSFEM. As a result, the corresponding stability conditions and the existence and 

uniqueness of the solution of the p x /i-SHSFEM scheme (4.26) follow from those of the 

k x /r-SHSFEM (cf. Lemma 4.2, Lemma 4.4 and Theorem 4.1). 

Following the same routine as in Section 4.2.3 (cf. the estimates (4.12)-(4.14)), we only 

need to estimate the terms inf ||<tv — and inf lu/v — v r ,Jr- Since 

T ph eZP®z Dh " ph "° Vph ez p®v Dh ' phh 


inf | |crjv - T ph \k 
T ph eZP®-E Dh 


< 


< 


inf 


+ 


inf 


I a N Tpno 

\CTN - Tp || 5 + h\\o- N \\j, 


t p £ZP®X d " u T h &Ll(T)®Y, Dh 

inf 

TpGZP(g)E D 


\<TN r h \\ Q 

(4.27) 


inf |UjV Vp/,|r < inf |uAT-V p |r + inf Iwy-V/Jr 
Vph&ZP®V Dh F v p £Z P®V D F v h £L 2 p (r)®V Dh 

~ c l n P L/ l U iV-VplT + /i||Uiv|l2, 

~Vp^Z^<S>VD 


it remains to estimate inf ||<tjv — "tv.i in 

t p £Zp®i: d y u v p ezp®v D 

<=lZ P n\ we easily have the following estimates: 


(4.28) 

and inf |ujy — v p |y. Recalling Z p = 


inf llcrjv — T, 

o€.Zp<S)Y!i]) 


N 

PM 0 ~y] 

“Tp„ezg»®ED 
N 


\ a N — 'rp„||c' 0 (r,E £ ,)i 
U N - V Pnllc°(r,VD)- 


(4.29) 


,J D ,U |UW - V -'l ~ E v „„ e §Lvb l|UN ' (4.30) 

Then the thing left is to estimate the right hand side terms of the above two inequalities. 

Denote T* := r *> then r = r n x T* , and for any y G T we denote y = (y n , y*) 

with y n G r n and y* G T*. We have the following lemma. 

Lemma 4.5. Let (cr/v, u jy) G (L^(r) < 8 > E^) x (L 2 p (T) < 8 > Vo) be the solution of the prob¬ 
lem (3.21). Then for any x G D, y = (y n , y* n ) G F, the solutions <T/v(x, y n ,y*) and 
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UN(x,yn,Vn) as functions of y n , i.e. a N : T n C^r^E^), u N : T n C°(T* n ;V D ), can 
be analytically extended to the complex plane 

dn) .— { z G C, dist(z, r n ) ^ d n f, 


with 0 < d n < 2 ~ and y n given by (4.17). In addition, for all z G E(T n ;d n ), it holds 
ll <J Af( 2: )llc 0 (r*;S) + l' u Jv( z )lc 0 (r*;VD) ~ T~oTir d l/zv lie 0 (t-,l 2 (d)) +11 ffiv 11 c°(r ; L 2 (aDi )) + !)• 

' i za n 7 n 

(4.31) 


Proof. Similar to (4.16), for y G f, r > 0 and n = 1,2, ...,N it holds 


liar o-jv(-,y)|| 0 , |ai n u JV (-,y)|i 


(27n) r (||fjv(-,y)||o + ||gjv(-.y)llo,ar» 1 + !)• (4-32) 


For any y n G F n , we define power series 

OO / \T* OO / \ T* 

^jv(x, 2 , y* n ) = ^ n <rjv(x, y n , y*n ), uat(x, z, y* n ) = a^ujv(x, y n , y* 


then it follows 


00 I 7 _ n, \r 

kiv(x,2,y*)||o 7 — l|ay„gjv(x,y n ,y*)||o, 


oo I I 

|ujv(x,Z,y*)|i |^„uiv(x,2/n,yn)|l- 

r=0 

Due to (4.32), we easily know that the above two series converge for all z G H(r n ; d n ). Fur¬ 
thermore, by a continuation argument, the functions er^y, ujv can be extended analytically 
on the whole region E(r n ;d n ), and the estimate (4.31) follows. □ 


In order to estimate the right-hand-side terms of (4.29) (4.30), we need one more lemma 
by Babuska et al [2], 


Lemma 4.6. Let B be a Banach space, and L C R be a bounded set. Given a function v G 
C°(L-,B) which admits an analytic extension in the region of the complex plane S(L;d) = 
{z G C, dist(z, L) < d } for some d > 0, it holds 

2 

min \\v-w\\ c o, L . B) < - -g~ p max ||v(z)||s, (4.33) 

w&P p (L)®B ' Q — 1 z£E(L-,d) 

where P p (L ) := span(y s , s = 0,1, 1 < g := ^ + y^l + j j^. 
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In light of (4.27)-(4.30) and Lemmas 4.5-4. 6 , we immediately obtain the following 
result. 

Theorem 4.4. Let (ctn,un) € {L 2 (T) <g> E&) x (L 2 ( r) ® Vp) and (cr p h> u pk ) G (Z p <g> 
x {Z p ®Voh) be the solutions of (3.21) and (4-26), respectively. Then, under the 
same condition as in Lemma 4-1 cmd for sufficiently large N, it holds 

N 

||crjv - CTkhWo + \u N ~ U kh Iy ^ h + y Qn~ Pri , (4.34) 

n =1 

where g n = ^ + y/l + and 0 < d n < 

The above theorem, together with Lemma 3.4, implies the following a priori error 
estimates for the p x h-SHSFEM approximation (cr p h, u p h). 

Theorem 4.5. Let(a,u ) G (Lp(f2; Ep))x(L 2 p (Ll-, Vo)) and (a pk ., u pk ) € {Z p ®Ep h , Z p ® 
Voh) be the solutions of (2-4) and (4-26), respectively. Then, under the same conditions 
as in Theorem 4-4> it holds 

N 

Ik - CkhW-Q + |m- Ufcfcli < A rl/2 e _rArl/2 + /1 + ^Qn~ Pn (4.35) 

n= 1 

/or any r > 0 if the covariance functions of E, f and g are piecewise analytic, and holds 

N 

Ik - a-kh\\o + \u- u kh | T < N~ s + h + ^2 2n~ Pn (4.36) 

n =1 

for any s > 0 if the covariance functions of E, f and g are piecewise smooth. 

Remark 4.3. This theorem shows the p x h-SHSFEM yields exponential rates of conver¬ 
gence with respect to the degrees (pi,P 2 , ...,Pn) of the polynomials used for approximation. 


5 Numerical examples 


In this section we compute two numerical examples to test the performance of the 
proposed p x li -version of stochastic hybrid stress finite element method. We note that the 
p x /r-SHSFEM can be viewed as a particular case of the k x h version. For convenience 
we denote 


0-u 


U - u/dy 



C(j . 


<y - cr h 



o 


22 




Rectangular meshes 



10 x 2 


Irregular meshes 


Figure 2: Finite element meshes 


where (u h-i&h) is the corresponding stochastic finite element approximation to the exact 
solution (u, <t). 

Example 1 : stochastic plane stress problem 

Set the spatial domain D = (0,10) x (—1,1) with meshes as in Figure 2. The body force 
f and the surface traction g on dD\ = {(xi,X 2 ) € [0,10] x [—1,1] : x\ = 10 or X 2 = ±1} 
are given by 

f=(0,0) T , gUi=io = (-2Ex 2 ,0) T , g|x 2 =±i = (0,0) T . 

The exact solution (u, a) is of the form 

—2xiX2 \ / —2Ex2 0 

x\ + v{x\ — i) J y o o 

where E is a uniform random variable on [500,1500], and we set v = 0.25. 

In the computation we use the exact form of the stochastic coefficient E and take 
N = 1, so there is no truncation error caused by the K-L expansion in the approximation. 
Numerical results at different meshes and different values of p are listed in Tables 1-2. 
For comparison we also list results computed by a stochastic finite element called PC x h 
method, where the polynomial chaos (PC) method [9] and the PS element method are used 
in the stochastic field T and the space domain D , respectively. In the PC x h method, 
p denotes the degree of polynomial chaos. We note that the computational costs of the 
PC x h method and the p x /r-SHSFEM are almost the same with the same p. 

From the numerical results we can see that the solutions are more accurate with the 
increasing of p and the refinement of meshes. Especially, p = 1 and p = 2 for the p x h- 
SHSFEM give almost the same results, which implies that the solutions are accurate 
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enough with respect to the p -version approximation of the stochastic field for given spatial 
meshes; In these cases, the px/t-SHSFEM is of first order accuracy in the mesh size h for the 
displacement approximation and yields quite accurate results for the stress approximation. 
What’s more, we can see that the p x h-SHSFEM is more accurate than the PC x h method 
at the same p. 


Table 1: 

Results for two methods under rectangular meshes: 

: Example 1 

Methods 



&IL 







P 

5x1 

10 x 2 

20 x 4 

40 x 8 

5x1 

10 X 2 

20 x 4 

40 x 8 


4 

0.0733 

0.0375 

0.0204 

0.0130 

0.0202 

0.0202 

0.0202 

0.0202 

PC xh 

6 

0.0728 

0.0365 

0.0186 

0.0098 

0.0079 

0.0079 

0.0079 

0.0079 


8 

0.0727 

0.0364 

0.0182 

0.0092 

0.0033 

0.0033 

0.0033 

0.0033 


0 

0.1223 

0.1050 

0.1003 

0.0990 

0.2774 

0.2774 

0.2774 

0.2774 

p x h 

1 

0.0727 

0.0363 

0.0182 

0.0091 

0 

0 

0 

0 


2 

0.0727 

0.0363 

0.0182 

0.0091 

0 

0 

0 

0 

Table 2 

: Results for two methods under irregular meshes: 

Example 1 




&IL 




^<7 



Methods 

P 

5x1 

10 x 2 

20 x 4 

40 x 8 

5x1 

10 x 2 

20 x 4 

40 x 8 


4 

0.1431 

0.0637 

0.0325 

0.0181 

0.2632 

0.0579 

0.0231 

0.0203 

PC xh 

6 

0.1429 

0.0631 

0.0314 

0.0160 

0.2626 

0.0549 

0.0137 

0.0083 


8 

0.1429 

0.0630 

0.0312 

0.0156 

0.2625 

0.0544 

0.0117 

0.0041 


0 

0.1435 

0.1160 

0.1037 

0.0999 

0.3684 

0.2816 

0.2775 

0.2774 

p x h 

1 

0.1429 

0.0630 

0.0311 

0.0155 

0.2524 

0.0509 

0.0104 

0.0023 


2 

0.1429 

0.0630 

0.0311 

0.0155 

0.2524 

0.0509 

0.0104 

0.0023 


Example 2 : stochastic plane strain problem 

The domain and meshes are the same as in Figure 2. The body force f = (0,0) T . 
The surface traction g on dD\ = {(aq, aq) £ [0,10] x [—1, 1] : aq = 10 or aq = ±1} is given 
by gU 1= io = (— 2Eaq, 0) t , g\ X2 =±i = (0,0) T , and the exact solution (u, a) is of the form 

— 2(1 — v 2 )x 1X2 \ ( —2Ex2 0 

(1 — v 2 )x\ + v{l + v)(x 2 — 1) J y 0 0 

where E = 1 + £ 2 , £ is a standard normal Gaussian random variable. 

Similar to Example 1, in the computation we use the exact form of the stochastic 
coefficient E and take N = 1. Numerical results at different meshes, different values of 
p and different values of Poisson ratio u are listed in Tables 3-8. For comparison we also 
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list results computed by a stochastic finite element called pxbilinear method, where the 
p-version method and the bilinear element are used in the stochastic field T and the space 
domain D, respectively. We note that the computational costs of the px bilinear method 
and the p x /i-SHSFEM are almost the same. 

Tables 3-4 show that the pxbilinear method deteriorates as v —> 0.5 or A —> +oo, 
while Tables 5-8 show that the p x /i-SHSFEM yields uniformly accurate results for the 
displacement and stress approximations. Moreover, p = 0 and p = 2 give almost the same 
results, which implies that the solutions are accurate enough with respect to the p-version 
approximation of the stochastic field for given spatial meshes. 


Table 3: Results of e u for Example 2: pxbilinear method, p = 0 




Rectangular 

meshes 



Irregular 

meshes 


V 

10 x2 

20x4 

40x8 

80x16 

10 x2 

20x4 

40x8 

80x16 

0.25 

0.5384 

0.3061 

0.1625 

0.0883 

0.6854 

0.4501 

0.2532 

0.1356 

0.49 

0.8516 

0.6523 

0.4034 

0.2175 

0.8782 

0.7424 

0.5218 

0.3038 

0.499 

0.9533 

0.9070 

0.7856 

0.5579 

0.9511 

0.9145 

0.8322 

0.6617 

0.4999 

0.9661 

0.9556 

0.9365 

0.8760 

0.9641 

0.9550 

0.9378 

0.8925 



Table 4: Results of 

e u for Example 2: 

pxbilinear method, 

p = 2 




Rectangular 

meshes 



Irregular 

meshes 


V 

10 x2 

20x4 

40x8 

80x16 

10 x2 

20x4 

40x8 

80x16 

0.25 

0.5384 

0.3061 

0.1625 

0.0883 

0.6854 

0.4501 

0.2532 

0.1356 

0.49 

0.8516 

0.6523 

0.4034 

0.2175 

0.9511 

0.9145 

0.8322 

0.6617 

0.499 

0.9533 

0.9070 

0.7856 

0.5579 

0.9511 

0.9145 

0.8322 

0.6617 

0.4999 

0.9661 

0.9556 

0.9365 

0.8760 

0.9641 

0.9550 

0.9378 

0.8925 



Table 5: Results of e u for Example 2: p x h SHSFEM, p 

= 0 




Rectangular 

meshes 


Irregular 

meshes 


V 

10 x2 

20x4 

40x8 

80x16 10x2 

20x4 

40x8 

80x16 

0.25 

0.0372 

0.0186 

0.0093 

0.0046 0.0676 

0.0323 

0.0158 

0.0079 

0.49 

0.0488 

0.0244 

0.0122 

0.0061 0.0763 

0.0371 

0.0183 

0.0091 

0.499 

0.0497 

0.0248 

0.0124 

0.0062 0.0770 

0.0375 

0.0185 

0.0092 

0.4999 

0.0497 

0.0249 

0.0124 

0.0062 0.0770 

0.0375 

0.0185 

0.0092 
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Table 6: Results of e a for Example 2: p x h SHSFEM, p = 0 




Rectangular 

meshes 



Irregular 

meshes 


V 

10x2 

20x4 

40x8 

80x16 

10x2 

20x4 

40x8 

80x16 

0.25 

0 

0 

0 

0 

0.1513 

0.0866 

0.0450 

0.0227 

0.49 

0 

0 

0 

0 

0.1559 

0.0877 

0.0451 

0.0227 

0.499 

0 

0 

0 

0 

0.1563 

0.0878 

0.0452 

0.0227 

0.4999 

0 

0 

0 

0 

0.1564 

0.0878 

0.0452 

0.0227 


Table 7: Results of e u for Example 2: p x h SHSFEM, p = 2 




Rectangular 

meshes 



Irregular 

meshes 


V 

10x2 

20x4 

40x8 

80x16 

10x2 

20x4 

40x8 

80x16 

0.25 

0.0372 

0.0186 

0.0093 

0.0046 

0.0676 

0.0323 

0.0158 

0.0079 

0.49 

0.0488 

0.0244 

0.0122 

0.0061 

0.0763 

0.0371 

0.0183 

0.0091 

0.49 

0.0497 

0.0248 

0.0124 

0.0062 

0.0770 

0.0375 

0.0185 

0.0092 

0.4999 

0.0497 

0.0249 

0.0124 

0.0062 

0.0770 

0.0375 

0.0185 

0.0092 


Table 8: Results of e a for Example 2: p x h SHSFEM, p = 2 




Rectangular 

meshes 



Irregular 

meshes 


V 

10x2 

20x4 

40x8 

80x16 

10x2 

20x4 

40x8 

80x16 

0.25 

0 

0 

0 

0 

0.1513 

0.0866 

0.0450 

0.0227 

0.49 

0 

0 

0 

0 

0.1559 

0.0877 

0.0451 

0.0227 

0.499 

0 

0 

0 

0 

0.0156 

0.0878 

0.0452 

0.0227 

0.4999 

0 

0 

0 

0 

0.1564 

0.0878 

0.0452 

0.0277 
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